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The objective of this article is to use Fourier-Hankel decomposition as suggested earlier by Braun & Fan (1998) to estimate 
the integrated horizontal meridional flow velocity as a function of mode penetration depth, and to find ways of potentially 
improve this technique. We use a time series of 43200 (30 days) consecutive full-disk Dopplergrams obtained by the MDI 
(Michelson Doppler Imager) instrument aboard the SOHO (Solar Heliospheric Observatory) spacecraft in April 1999. 
We find averaged meridional flow estimates of 15 m/s for modes with a penetration depth in the upper 20 Mm of the 
solar convection zone. This reproduces the results of the earlier investigations. Moreover we conclude that this method 
has the potential to become a new technique to measure the meridional circulation in the deep convection zone, if some 
improvements will be applied. 
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On the solar surface a meridional flow of 10-20 m/s can be 
measured. This flow is directed away from the equator to- 
wards the poles (Hathaway (1996), Komm, Howard & Har- 
vey (1993),Latushko (1996), Snodgrass & Daily (1996)). 
Because of mass conservation a poleward meridional flow 
in the outer part of the convection zone requires an equa- 
torward return flow in deeper layers. This meridional cir- 
culation might play an important role in the dynamo pro- 
cess which causes the solar magnetic cycle (Choudhuri, 
Schussler & Dikpati (1995),Dikpati & Gilman (2006),Dik- 
pati & Charbonneau (1999), Brandenburg, Moss & Tuomi- 
nen (1992)). It has been the topic of various studies (starting 
with Giles et al. (1997)) to infer the depth-dependent profile 
of this meridional circulation with various techniques from 
helioseismology. 

The motivation for this work is to reproduce a method 
suggested by Braun & Fan (1998) for inferring the speed of 
the averaged solar meridional flow. The method is based on 
splitting up the oscillation signal observed in a time series 
of Dopplergrams into two wave-fields of anti-parallel travel- 
ling waves by a Fourier-Hankel decomposition (Braun, Du- 
vall, & Labonte (1987)). Originally an average meridional 
flow in the order of 10 m/s was found in the upper part of 
the convection zone. We are able to reproduce the results 
on an independent data set and with a newly created data 
analysis pipeline. The reproduction of the method provides 
a first step in measuring the meridional circulation. We dis- 
cuss potential improvements of this method to obtain more 
accurate results and to possibly reach to greater depths in 
future studies. 



2.1 Normal mode decomposition 

Braun, Duvall, & Labonte (1987) applied Fourier-Hankel 
spectral analysis to decompose locally the solar oscilla- 
tion signal into wave fields travelling in opposite directions. 
An overview to the technique is given in Gizon & Birch 
(2005). This decomposition was used for investigating the 
influence of sunspots on the helioseismic p-modes. It al- 
lowed detecting phase-shifts and damping of p-modes by 
sunspots (Braun, Duvall & Labonte (1988), Bogdan et al. 
(1998), Braun (1995)). 

A further application allowed studying the effect of the 
meridional flow on the solar oscillations (Braun & Fan 
(1998)). However, the Fourier-Hankel decomposition was 
originally designed for studying signals on a small section 
of a spherical surface by approximating the surface by a pla- 
nar geometry. Therefore it is generally not suited for study- 
ing global flows on a sphere. Nevertheless, we reproduce the 
data analysis in the following. But we note, this method can 
only be seen as an initial approximative step to obtain rough 
estimates of average effects. It might allow to conclude on 
what to expect from a more proper analysis in future. 

The velocity amplitude signal of the solar eigenoscil- 
lations ^(6,ip,t), consisting of standing p-modes, can be 
locally separated into two fields of travelling waves on the 
sphere. This local decomposition is carried out within an 
annular region around a central point of interest. This point 
is the piercing point of the axis of the used spherical coordi- 
nate system, where 9 is the latitude, and <p is the longitude. 
In the original concept (Braun, Duvall, & Labonte (1987)) 
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this point has been the center of a sunspot or an active re- 
gion. In the following we investigate a subsection of an an- 
nular region around the poles of the sun, i.e. the polar axis 
of the coordinate system is identical with the rotation axis 
of the sun. The coordinates used therefore coincide with the 
standard spherical coordinate system. 

According to Braun & Fan (1998) the oscillation veloc- 
ity signal W is expanded as 

*(9,<p,t) = e i(mv+27ri/t) (1) 

I m is 



x [Ai mv er (cos 9) + B lmu (er)* (cos 9)} 



with 

er(cosfl) :=N? 



2i 



Pr (cos 9) + -QT (cos 9) , (2) 

IT 

where the asterisk denotes complex conjugation, P™ and 
Q™ are Legendre functions of the first and second kind, 
and iVj m is a normalisation constant. The complex quanti- 
ties Ai mu and B\ mv are the amplitudes of the inward and 
outward going waves respectively, n is the radial order, I is 
the harmonic degree, m is the azimuthal order, t is time and 
v is the temporal frequency. We refer the reader to Braun & 
Fan (1998) and Braun, Duvall & Labonte (1988) for more 
details on these decomposition. 

Following Braun & Fan (1998) Hankel functions are 
used as approximations to (P™ ± 2i/7rQ™) 
Xl-m)\ 



H^(L9) « (-iy 



(l + m)\ 



2i 



P{ n (cos9) ± -Qf (cos 9) 



(3) 



with L — yjl(l + 1). This approximation is valid in the 
limit / 3> m, which is fulfilled in the current investiga- 
tion (Braun (1995)). A decomposition into wavefields by 
applying a method built on Legendre functions, rather than 
Hankel functions, would be more accurate for the spherical 
case. But the application of Hankel functions provides an 
estimation which is much faster to calculate. 

According to Gizon & Birch (2005) in some cases a dis- 
crete set of Hankel functions can be selected which guaran- 
tees orthogonality. Explicitly an integration over the entire 
range of ip would provide such an orthogonal set of func- 
tions in the present case. In the further data analysis, how- 
ever, the integral is carried out only over a limited range 
in ip. This then introduces possibly leakage effects between 
modes of different azimuthal order m. Therefore only the 
best possible approximation to an orthogonal set can be ac- 
tually selected. 

The wave amplitudes Ai mv and Bi mv are extracted from 
the wave field ^(9, ip, t) according to Braun et al. (1992) by 



A 



C 



^(9,ip,t)H^ 2 \L j 9) 



xe -(^+^t) $ded(pdt ) 
Bi mv ~cJ ^>(9,<p,t)H^(L j 9) 



(4) 



(5) 



where C is a normalisation which is given approximately by 
L/ (4T0) with the range in latitude over which the inte- 
gration is carried out, and T the duration of the observation. 

Summarized, the signal is decomposed by a spatial Han- 
kel transformation and a temporal Fourier transformation. 
The result of this decomposition is a set of frequency spec- 
tra \Ai mL ,\ 2 and |P; m „| 2 depending on the harmonic degree 
I. Peaks occurring in these spectra are to be identified as 
modes of different radial order n with the frequencies 
and iA respectively. 

Potentially differences between the mode frequencies 
for poleward and equatorward propagating waves with same 
indices emerge due to advection. The frequency differences 
carry information about the influence of the horizontal com- 
ponent of the meridional circulation on the p-modes. Thus 
an exploration of the set of spectra leads to a set of fre- 



quency shifts Av n i 
of unique indices. 



assigned to a respective set 



2.2 Effect of meridional circulation on the p-modes 

In the following we provide a formula that describes the 
effect of a horizontal flow U on the p-mode frequencies. 
For detailed steps to obtain this formula see Braun & Fan 
(1998) and references therein. 

Under the influence of a horizontal flow U the net 
frequency shift between waves propagating poleward and 
equatorward is given by 



/ / ((U) e /r)K nl (r) dr 
Au nl = -2— 

Hq 

tt / K n i(r) dr 



(6) 



where K n i (r) is a depth- and mode-dependent weighting 
function and 



(U)e ■= 



"max 



dO. 



(7) 



is the flow average over the spanned polar angle. Follow- 
ing Braun & Fan (1998) the frequency shift of a zonal mode 
due to advection is proportional to the horizontal wavenum- 
ber and the quantity 

(U) := Au ln nR Q /l , (8) 

gives the weighted depth average of the angular mean 
meridional flow (U)e/r multiplied by the solar radius. 
Together with 
1 



v :- 



o Kl +V nl) . 



(9) 



all these considerations yield a data set {(({/), 
where each i corresponds to a given p-mode. 

In the last step the turning points of the respective modes 
are determined by using the relation 

c 2 (r t ) (2TTP) 2 



1(1 + 1) 



(10) 



© 0000 WILEY- VCH Verlag GmbH & Co. KGaA, Weinhdm 



www.an-journal.org 



Astron. Nachr. / AN (0000) 



3 



with r t the depth of the inner turning point. From this we 
obtain the mode penetration depth e := R@ — r t below the 
solar surface. 

We note that the estimation of the integral carried out is 
valid only in very shallow layers as it required approximat- 
ing (U)e/r by (U)$/Rq. The velocity amplitude inferred 
later is therefore only reliable within the outer layers of the 
Sun. In deeper layers, however, this estimation might only 
provide information about the quality of the velocity profile, 
i.e. the sign of the weighted velocity average. 

2.3 Data analysis 

The analysed data consist of a time series of full disk Dopp- 
lergrams recorded by SOHO/MDI. The timeseries lasts 
from April 1-30, 1999 and consists of 43200 consecu- 
tive single Dopplergrams. Each Dopplergram is given on a 
1024x 1024 pixel grid. The duty cycle of the time series is 
better than 99%. The data used by Braun & Fan (1998) were 
timeseries of MDI and GONG Dopplergrams from 1997. 

Each Dopplergram is projected onto an equidistant d-ip- 
lattice, referring to the standard three-dimensional spherical 
coordinate system, ^-centered at the solar north pole. After 
the transformation the Dopplergram is asymmetric with re- 
spect to the equator, because at the observing time the north- 
ern hemisphere of the Sun was better visible. The oscillation 
signal is therefore better measurable on the northern hemi- 
sphere. 

Two areas on the solar surface are investigated, one in 
the northern the other in the southern hemisphere. The two 
investigated zones are placed symmetrically to the equa- 
tor with a polar interval of = ir/4 and an azimuthal 
interval of ipm = tt/2. The selected areas are centered at 
{6> = tt/4,</3 = tt/2} and {9 = 3/4 n,(p = tt/2} re- 
spectively. The sizes of the fields were chosen such to yield 
the best oscillation signal-to-noise ratio. Moreover this al- 
lows the resolution of oscillations with low harmonic de- 
gree I in order to reach to greater depths in the solar interior. 
This choice is comparable with the investigation carried out 
by Braun & Fan (1998). 

The geometry of the chosen fields restricts the possi- 
ble resolution of the harmonic degree to I > 4 as a lower 
boundary. The resolution of the Dopplergrams yields an up- 
per boundary of I < 1024. As we assume a (,5-independent 
meridional circulation, only modes with m = are investi- 
gated. 

3 Results 

3.1 General results for both hemispheres 

The result of our investigation are measurements of the av- 
erage horizontal meridional velocity (U) as a function of 
v/L. Besides the relation (U) ~ v/Lwe are interested in 
the average horizontal velocity as a function of the penetra- 
tion depth e. FigureQ]displays the general results. 
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Fig. 1 Estimated averaged horizontal meridional velocity 
on the northern (top) and southern (bottom) hemisphere as a 
function of v/L and as function of mode penetration depth 
e. 



Because of the large scatter in the results we applied a 
binning of the resulting data points. In this particular case, a 
binning window of five points has been chosen. This choice 
reflects a trade-off between a desired rather smooth profile 
on the one hand and the resolution, which on the other hand 
should not be too coarse. The reported error bars given in 
the plots below are due to the statistical scatter of the points 
within one binning-window. In the plots a positive value of 
the velocity is equivalent to a poleward flow on the respec- 
tive hemisphere. 

In particular, Fig. Q] (top) gives the averaged horizon- 
tal meridional flow on the northern hemisphere. Obviously, 
various systematic trends are noticeable. First, for modes 
penetrating between 1-4 Mm the results are distributed 
around m/s. This feature has not been seen in other in- 
vestigations of the meridional flow and may correspond to 
effects of granulation and supergranulation. The true ori- 
gin of this is, however, not clear. Second, for modes prob- 
ing deeper below the solar surface there is a clear positive- 
valued trend of (15±7) m/s. The error bars in this region are 
relatively small. Third, the observed variation in the curve 
might be an artefact due to different sensitivities of the p- 
modes to the flow. For the given mode penetration depths, 
the results agree well with meridional flow measurements 
by other local helioseismology techniques like ring-diagram 
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analysis (Haber et al. (2002), Zaatri et al. (2006)) or time- 
distance helioseismology (Zhao & Kosovichev (2004)). 

Figure Q] (bottom) displays the results for the southern 
hemisphere. As in the former case, for modes turning close 
to the solar surface the points are distributed around m/s. 
For modes turning further down in the deep interior a similar 
general trend of the curve is observable as on the northern 
hemisphere. However, the average velocity is about m/s, 
in contrast to the northern hemisphere where a value of ca. 
(15 ± 7) m/s is obtained. Additionally the curve shows a 
far larger amount of scatter, not only expressed in terms of 
larger error-bars but also in the shape of the curve itself. One 
reason for this difference between the northern and south- 
ern hemisphere might be an error in the determination of 
the direction of the rotation axis of the Sun (Giles (2000)). 
Another reason might arise by leakage effects from modes 
with m 7^ whose frequencies are shifted by differential 
rotation. In the following steps we further analyse the mea- 
surements for the northern hemisphere data only. 

We compare our results obtained for the northern hemi- 
sphere with the results presented by Braun & Fan (1998). 
We analyse a data set set from 1999 whereas Braun & Fan 
(1998) analysed a data set from 1997. Both results agree 
well. Differences can only be found on small scales. We 
considered only modes with an azimuthal order of m = 0, 
whereas the comparable work used a bigger set with m/0 
and applied an additional averaging. Therefore our errorbars 
are bigger. Overall, this gives confidence in our newly de- 
veloped data analysis technique and the resulting measure- 
ments. 

3.2 Possible extension to greater depths 

We are able to measure frequency shifts for modes with 
harmonic degrees of I = 73 - 1013. The modes with the 
lower harmonic degree penetrate to depths of about 150 
Mm. Therefore this method might allow to derive informa- 
tion about the meridional flow throughout 75% of the con- 
vection zone. Fig. [2] gives the averaged horizontal merid- 
ional velocity (U) estimated by this method for modes pen- 
etrating to such depths. Between the values of e = 20 Mm 
and e = 100 Mm there is no clear behaviour visible. Be- 
low e ~ 110 Mm the curve shows a trend to become 
negative-valued. However there are only a few measure- 
ments at great depths. Moreover the results are largely scat- 
tered which leads to the respective large error-bars. Never- 
theless, the measured frequency shifts contain information 
on the meridional flow down to approximately 150 Mm. 
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Fig. 2 Estimated integrated horizontal meridional veloc- 
ity on the northern hemisphere of the Sun as a function of 
v/L and mode penetration depth e. 

Based on our results we are able to suggest improve- 
ments for this technique which we discuss in the following 
section. 

4 Discussion and Future Work 

We applied a Fourier-Hankel decomposition to MDI Dopp- 
lergrams as described by Braun & Fan (1998) to estimate 
the average effect of the meridional flow on solar p-modes. 
We find an average poleward meridional flow velocity on 
the northern hemisphere in the order of 15 m/s which af- 
fects modes penetrating in a depth range of 4-20 Mm. Oc- 
curring advection effects on modes penetrating to shallow 
layers (e « 1-4 Mm) are interpreted as the effect of (super- 
) granulation. However the latter feature was not found in 
other helioseismic investigations; its true origin needs to be 
clarified. Furthermore frequency shifts for modes with an 
inner turning point at about 150 Mm were measured. 

It is encouraging that for modes turning in the up- 
per 20 Mm of the convection zone the inferred average 
meridional flow agrees well with other independent mea- 
surements of the meridional flow velocity (Braun & Fan 
(1998), Giles (2000), Haber et al. (2002), Zhao & Koso- 
vichev (2004), Zaatri et al. (2006)). Furthermore, as there 
is an almost perfect agreement between our results and the 
results of Braun & Fan (1998) we are sure that we were able 
to rebuilt their original concept. 

The clear differences between the northern and southern 
hemispheres could be either due to unresolved asymmetries 
in the observational setup at the given observing time or due 
to other systematic errors. At least for MDI data an offset of 
typically about 10 m/s in measurements of the meridional 
flow between the two hemispheres occur, independent from 
the method of investigation (see e.g. Braun & Fan (1998) or 
Zaatri et al. (2006)). 

The presented method is only a rough first-step- 
estimation of the profile of the meridional circulation. Fu- 
ture investigations using the decomposition of the wavefield 
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in space and time and aiming at higher accuracy in the re- 
sult should take into account the spherical symmetry of the 
problem. Using Legendre functions instead of Hankel func- 
tions will be a first but vital step. Additionally the effect of 
leakage, potentially not only based on a finite set of degrees 
(c.f. Gizon & Birch (2005)) but also on the neglected influ- 
ence of modes with m/0, has to be analysed. 

We conclude that the presented method - as is - pro- 
vides a tool for an easy-to-handle analysis of the average 
solar meridional flow profile in shallow layers. Using longer 
time series, this method might allow a better frequency res- 
olution to obtain better estimates for the frequency shifts. 
Moreover, there is some freedom in choosing the investi- 
gated range over polar angle. This might allow to assess 
estimates for the average meridional flow as a function of 
latitude, too. 

The used assumptions break down for modes penetrat- 
ing to layers below 20 Mm. Moreover the scatter in the data 
is high at great penetration depths. Nevertheless, we find a 
negative-value trend in the estimation of the average merid- 
ional flow for p-modes probing deeper than 110 Mm. Such 
a negative value might give a hint to an equatorwards di- 
rected meridional flow in the layers above. However, before 
concluding on the meridional return flow all entrapments in 
the data analysis need to be resolved. 

Most important, a proper inversion of the measured fre- 
quency shifts needs to be carried out to estimate the veloc- 
ity profile over a certain depth range. This will require to 
deal with the forward problem properly. A normal mode ap- 
proach needs to be carried out to obtain not only a proper es- 
timation of the effect of all three spatial components of the 
meridional flow on the mode frequencies but also to obtain 
integral kernels. Then it will be interesting to see, whether 
enough data and information could be collected to construct 
averaging kernels at greater depths. If all these suggested 
improvements are incorporated, then the meridional flow 
could be studied in better detail as function of depth, lati- 
tude and time with this method. 
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